Tuning Jammed Frictionless Disk Packings from Isostatic to Hyperstatic 
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We perform extensive computational studies of two-dimensional static bidisperse disk packings 
using two distinct packing-generation protocols. The first involves thermally quenching equilibrated 
liquid configurations to zero temperature over a range of thermal quench rates r and initial packing 
fractions followed by compression and decompression in small steps to reach packing fractions 4>j at 
jamming onset. For the second, we seed the system with initial configurations that promote micro- 
and macrophase-separated packings followed by compression and decompression to <f>j. Using these 
protocols, we generate more than 10 4 static packings over a wide range of packing fraction, contact 
number, and compositional and positional order. We find that amorphous, isostatic packings exist 
over a finite range of packing fractions from </> m i n < 4>J < </>max in the large-system limit, with 
0max ~ 0.853. In agreement with previous calculations, we obtain <p m i n « 0.84 for r > r* , where r* 
is the rate above which cf>j is insensitive to rate. We further compare the structural and mechanical 
properties of isostatic versus hyperstatic packings. The structural characterizations include the 
contact number, bond orientational order, and mixing ratios of the large and small particles. We 
find that the isostatic packings are positionally and compositionally disordered, whereas bond- 
orientational and compositional order increase with contact number for hyperstatic packings. In 
addition, we calculate the static shear modulus and normal mode frequencies of the static packings 
to understand the extent to which the mechanical properties of amorphous, isostatic packings are 
different from partially ordered packings. We find that the mechanical properties of the packings 
change continuously as the contact number increases from isostatic to hyperstatic. 

PACS numbers: 83.80.Fg61.43.-j62.20.-x64.75.Gh 



I. INTRODUCTION 



The ability to enumerate and classify all of the me- 
chanically stable (MS) packings of frictionless particles 
is important for understanding glass transitions l] in 
atomic, molecular, and colloidal systems, and the struc- 
tural and mechanical properties of particulate materials 
such as granular media, foams, and emulsions. For exam- 
ple, if all MS packings in a given system are known, one 
can measure accurately the frequency with which each 
MS packing occurs, and determine how the packing fre- 
quencies and materials properties depend on the prepa- 
ration history 0, |Jf. Further, MS packing frequencies 
are important for identifying the appropriate statistical 
mechanical ensemble for weakly perturbed granular ma- 
terials [<|. However, since the number of MS packings 
grows exponentially with the number of particles [5[ , ex- 
act enumeration of static packings is prohibitive for even 
modest system sizes [6]. Thus, one of the most important 
outstanding questions in the area of disordered particu- 
late materials is determining how the packing-generation 
protocol influences the distribution of MS packings and 
their structural and mechanical properties. 

Previous work has suggested that the positional order 
of MS packings of frictionless spheres increases mono- 
tonically with packing fraction and contact number in 
dense packings 0, Q. However, the MS packings in these 
previous studies were created using monodisperse sys- 
tems, which are prone to crystallization [9], and pre- 
pared using the Lubachevsky-Stillinger compression al- 



gorithm [lfj, which is a thermalized packing-generation 
protocol. In addition, these prior studies did not dis- 
tinguish the distribution of isostatic MS packings (in 
which the number of deg rees of freedom matches the 
number of constraints [ll|) from the distribution of hy- 
perstatic packings (with more contacts than degrees of 
freedom). Later work characterized bidisperse systems, 
which are less prone to crystallization, but focused on 
microphase-separated states, not amorphous, isostatic 
packings [12|]- However, recent studies on systems com- 
posed of 3D monodisperse, frictionless, spherical parti- 
cles have pointed out that amorphous, isostatic packings 
can exist over a finite range of packing fraction in the 
large-system limit, with no correlation between positional 
order and packing fraction 113, |k4J. Moreover, simula- 
tions [15| and experiments [16[ on two-dimensional sys- 
tems also suggest a finite range of jamming onsets rather 
than a single packing fraction in the large system limit. 

Further, the body of work on jammed particulate sys- 
tems has emphasized the concept of point J, i.e. that 
there is a single packing fraction at which jamming oc- 
curs in the large system limit [17l.ll8|. Since amorphous, 
isostatic packings can exist over a finite range of packing 
fractions, the onset of jamming should not be classified 
as a point in the jamming phase diagram, but rather as 
a region of finite extent. It has also been argued that 
the wide distribution of packing fractions at which the 
onset of jamming occurs in small periodic systems [17J ] is 
related to the finite range of packing fractions over which 
amorphous, isostatic packings occur in the large system 
limit [jjj. However, it has not been proved that these 
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two effects are directly connected. 

A number of overarching questions related to the con- 
nection between positional order, isostaticity, and mate- 
rial properties of static packings remain open. For ex- 
ample, can isostatic or nearly isostatic packings possess 
significant positional order and if so, what are the fun- 
damental differences in the normal modes and mechan- 
ical properties between those that do and do not pos- 
sess significant positional order? This question is par- 
ticularly important since recent studies have emphasized 
that amorphous, isostatic packin gs p ossess an excess of 
low-frequency normal modes [2; 
monic, ordered solids. 

In addition, previous work has drawn a strong con- 
trast between amorphous packings and configurations 
with crystalline order (22[. However, how different are 
the structural and mechanical properties of amorphous 
versus partially ordered particulate systems? For exam- 
ple, it is possible that the amorphous regions in the in- 
terstices between ordered domains in partially crystalline 
materials dominate the structural and mechanical prop- 
erties, in which case their properties would be similar to 
amorphous packings. At the very least, one would as- 
sume that there is not a strong difference between the 
mechanical properties of isostatic and only slightly hy- 
perstatic packings that possess significant positional or- 
der. 

In this article, we describe extensive computer simu- 
lations of collections of frictionless, bidisperse disks with 
short-range repulsive interactions to address two impor- 
tant, open questions: 1. What is the range of packing 
fractions over which amorphous, isostatic static packings 
occur with similar structural and mechanical properties, 
and 2. How do the structural and mechanical properties 
of static packings change with the deviation in the con- 
tact number at jamming onset from the isostatic value, 
zj — z iso [23]? Using two distinct packing-generation pro- 
tocols, we construct scatter plots for more than 10 4 static 
packings characterized by the contact number, packing 
fraction, measures of positional order, and mechanical 
properties. The first protocol involves thermally quench- 
ing equilibrated liquid configurations to zero temperature 
over a range of thermal quench rates r followed by com- 
pression and decompression in small steps to reach pack- 
ing fractions cj>j at jamming onset. For the second, we 
seed the system with initial configurations that promote 
micro- and macrophase-separated packings followed by 
compression and decompression to (j)j. 

Our main results are fourfold: 1. Isostatic, amorphous 
packings exist over a finite range of packing fraction from 
0min to i^max in the large system limit, with similar struc- 
tural and mechanical properties. 2. In agreement with 
previous calculations, we obtain </> m i n ~ 0.84 for r > r* , 
where r* is the rate above which </>j is insensitive to rate. 
In contrast, max depends sensitively on quench rate, sys- 
tem size, and boundary conditions. 3) The amorphous, 
isostatic packings coexist with an abundance of hyper- 
static, microphase- and macrophase-separated packings. 



4) When considering the full ensemble of static friction- 
less packings, the packings possess structural and me- 
chanical properties that span a continuous range from 
amorphous to partially ordered to ordered in contrast to 
the results and interpretations of recent studies [24j, [25[ . 
The remainder of the manuscript will be organized as 
follows. In Sec. mi we describe the computational system 
we consider and the two protocols we employ to generate 
static frictionless disk packings. In Sec. IIII1 we present 
our results, which include characterizations of the struc- 
tural (packing fraction, contact number, and several or- 
der parameters to detect positional and compositional 
order) and mechanical (shear modulus and eigenvalues 
of the dynamical matrix [3j) properties of more than 
10 4 static packings and comparisons of these properties 
for isostatic and hyperstatic configurations. Finally, in 
Sec. lIVi we provide our conclusions and promising future 
research directions. 



II. PACKING-GENERATION PROTOCOLS 

We focus on well-characterized two-dimensional sys- 
tems composed of N bidisperse disks (50-50 by num- 
ber), each of mass m, with diameter ratio d — o~i/o~ s — 
1.4 pjl, [It], [26| , within square, periodic simulation cells 
with side length L. We consider frictionless particles that 
interact through the finite-range, purely repulsive spring 
potential. The total potential energy per particle is given 
by 
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(1) 



where r^ is the center-to-center separation between disks 
i and j, e is the characteristic energy scale of the interac- 
tion, Q(x) is the Heaviside function, and oy = {<Ji+o-j)/2 
is the average diameter. We simulated a range of system 
sizes from N = 256 to 8192 particles to assess finite size 
effects. Energy, length, and time scales are measured in 
units of e, a s , and a s y/m/e, respectively. 

The packing fraction <pj at which jamming occurs and 
the structural and mechanical properties of static pack- 
ings can depend strongly on the packing-generation pro- 
tocol employed. Our goal is to generate static frictionless 
MS packings that span the range of contact numbers from 
the isostatic value Zi SO = 4 to the hexagonal crystal value 
-Zxtai = 6 and the range of positional order from amor- 
phous to phase-separated and from partially crystalline 
to crystalline states. To accomplish this, we investigate 
two distinct classes of packing-generation protocols: 1) 
thermal quenching from liquid initial conditions coupled 
with compression and decompression steps, which typ- 
ically generates amorphous configurations and 2) com- 
pression and decompression steps from initial conditions 
that promote micro- or macrophase separation 27]. 

Protocol 1: Thermal quenching from liquid initial con- 
ditions In this algorithm, we prepare equilibrated, liquid 



configurations at high temperature To = 10 -3 and in 
molecular dynamics (MD) simulations quench them to a 
very low final temperature Tf = 10~ 16 ~ at fixed pack- 
ing fraction 0.8 < fa < </>xtal — 7r/2\/3 [28| over a time 
interval t by rescaling the particle velocities so that the 
kinetic temperature T = TV -1 J^. mv 2 /2 obeys 



T(t) = T e- 



(2) 



where r is the thermal quench rate, which is varied over 
five orders of magnitude 10~ 5 < r < 1. We generated 50 
equilibrated, independent liquid configurations at T at 
each <pi by writing out configurations every lOr, where r 
is a decay time obtained from the self-intermediate scat- 
tering function at wavenumbers corresponding to the first 
peak in the structure factor .29] . 

After reaching a local potential energy minimum at 
each initial packing fraction fa and thermal quench rate 
r, we input the configurations into an 'athermal' algo- 
rithm ('packing finder') that searches for the nearest 
static packing in configuration space with infinitesimal 
particle overlaps. The algorithm has been described in 
detail in previous work [3|. Briefly, we successively in- 
crease or decrease the diameters of the grains (while 
maintaining the diameter ratio d), with each compres- 
sion or decompression step followed by conjugate gra- 
dient minimization of V. The system is decompressed 
when the total potential energy per particle at a local 
minimum is nonzero, i.e. there are finite particle over- 
laps. If the potential energy of the system is zero and 
gaps exist between particles, the system is compressed. 
The increment by which the packing fraction is changed 
at each compression or decompression step is gradually 
decreased. Numerical details of the algorithm are the 
same as in Ref. [3[. When this algorithm terminates, 
we obtain a static packing defined by the particle po- 
sitions {fi, r*2, . . . , rjv} and packing fraction <pj. Since 
we use an energy tolerance (per particle) V to i/e = 10~ 16 
for the termination of the energy minimization and com- 
pression/decompression scheme in the packing finder, the 
positions and packing fraction at jamming are extremely 
accurate with errors at one part in 10 8 . 

Protocol 2: Compression and decompression steps from 
initial conditions that promote order We will see below 
in Sec. [TTT] that Protocol 1 produces amorphous, isostatic 
packings. Thus, we seek an algorithm that will gener- 
ate static packings with variable positional and compo- 
sitional order. To bias the system toward micro- and 
macrophase-separated configurations, we seed the pack- 
ing finder with particular sets of initial conditions. We 
first divided the unit cell into s x s equal-sized parti- 
tions, where s is an even integer that ranged from 2 to 26, 
and placed approximately N / s 2 large or small particles in 
alternating partitions to create a checkerboard-like pat- 
tern. The particles were placed randomly in each parti- 
tion. The initial configuration is then input into the pack- 
ing finder to yield a static packing. In the large s limit, we 
expect amorphous static packings, while at intermediate 
and small s, we expect micro- and macrophase-separated 




FIG. 1: Average packing fraction (4>j) obtained from Proto- 
col 1 as a function of the negative logarithm of the thermal 
quench rate r for N = 1024. Data points at each rate repre- 
sent an average over typically 300 static, amorphous packings. 
The dashed line shows the scaling (cj>j) ~ [log 10 (r — ^*)] M , 
where /i ~ 0.5 and r* (a 0.03 is the thermal quench rate 
above which {(j>j) ~ 0.841 is independent of r. 



packings. To generate static packings near x tal we also 
divided the unit cell into two partitions and placed the 
large (small) particles on a hexagonal lattice in a region 
with area Al = d 2 /(l + d 2 ) (1 — Al) and then applied 
the packing finder. 



III. STRUCTURAL AND MECHANICAL 
PROPERTIES 

After generating static packings using the two packing- 
generation protocols described above, we contrast them 
by calculating several structural and mechanical proper- 
ties. The structural characterizations include the pack- 
ing fraction, contact number, and compositional and po- 
sitional order parameters. For the packing fraction at 
jamming onset, we calculate 
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1 



(3) 



including all N particles. For the contact number at 
jamming, we sum up all overlapping pairs (ry < o~ij) of 
particles, z.j = N c /N', where N' = N - N r , N r is the 
number of rattler particles with fewer than three con- 
tacts, and N c only includes overlapping pairs among the 
N' particles within the 'true' contact network. It is cru- 
cial to perform an error analysis on the contact number 
Zj, which is described in Appendix [Al 

Packing Fraction We show results for the average 
packing fraction (4>j) versus thermal quench rate r over 
five orders of magnitude obtained from Protocol 1 in 
Fig. Q] For large rates r > r* sa 0.03, the av- 
erage packing fraction (</>j) — > 0.841 is independent 
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FIG. 2: Scatter plot of the contact number zj versus the 
packing fraction at jamming onset <f)j. The open circles indi- 
cate static packings that were generated using Protocol 1 for 
N — 1024, while all other symbols indicate static packings 
generated using Protocol 2. The open squares, diamonds, 
and triangles correspond to N — 1024, 2048, and 4096, re- 
spectively, for all partitions s and systems with two parti- 
tions and random particle placements. The filled squares, dia- 
monds, upward triangles, and downward triangles correspond 
to N = 1024, 2048, 4096, and 8192, respectively, for the sys- 
tems with two partitions and initial crystal lattice positions. 
The black cross indicates the values zj = 6 and (f)j = 7r/2\/3 
for the hexagonal crystal. The labels (a)-(d) correspond to 
the images in Fig. [3] The inset shows the system-size de- 
pendence for systems with two partitions and random initial 
positions at N — 256 (leftward triangles), 1024 (squares), and 
4096 (upward triangles). 



of rate, which agrees with studies that employ at her- 
nial compression/decompression packing-generation al- 
gorithms |2j, |l7[. For r < r* , ((f) j) increases approx- 
imately as [log 10 (r — r*)] 05 with decreasing rate. We 
emphasize that all packings used to present the data in 
Fig. [1] are amorphous and isostatic. Since (<fij) increases 
so slowly, it is not possible to approach (/> x tai using pro- 
tocol 1. Using an extrapolation, we estimate that rates 
below 10 -45 are required to reach </> x t a i, and thus we em- 
ployed Protocol 2, not 1, to generate compositionally and 
positionally ordered packings. 

Contact Number In Fig. [5J we display a scatter plot 
of the contact number zj versus <pj for all static packings 
(where the contact number is insensitive to the definition 
of 'contact') generated using Protocols 1 and 2. (See Ap- 
pendix |A] for a discussion of the sensitivity of the contact 
number on the definition of contacting particles.) Fig. [5] 
shows several compelling features. First, nearly all of the 
static packings obtained from Protocol 1 (open circles) 
are isostatic with zj = 4, but they occur over a range of 
packing fractions m i n < 4>J < 0max, where (f> m i n = 0.837 
and </> max = 0.853. As shown in Appendix [X] max is 





"V^/V- 














JO'S 


iw^^SSfiffiC^f? 


Ti>V 








- . ■ 


< ^-ijerw^XhH i 


(c) 


sP^^^xy 








WMmt 




^^J^Vyx^cF^^ 


^S^kFS 




lllllll 




w8§fw 


-■ r03Q& 


wWlvSm 






r?i&fo$$$>n% 




(e) 






FIG. 3: Images of representative static packings from the 
scatter plot in Fig.[2]with (a) </>j = 0.837, zj = 3.99, (b) <f)j = 
0.853, zj = 4.00, (c) cj>j = 0.846, zj = 4.04, (d) cj>j = 0.860, 
zj = 4.41, and (e) <j)j = 0.892, zj ~ 4.1. (See Appendix^]) 
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FIG. 4: Scatter plot of the fraction of contacts between two 
large fu or two small particles f sa versus packing fraction <j>j 
for all static packings from both protocols. The diamonds 
(circles) and triangles (squares) display data from Protocol 1 
(2) for fu and f ss , respectively. 



likely only a lower bound for the largest packing fraction 
at which isostatic packings can occur in these systems. 
Second, we find a cluster of data points for Protocol 2, 
for which the average zj is strongly correlated — varying 
roughly linearly — with <pj. The cluster originates near 
4>j « 0.84, zj = Zi SO = 4. In the inset to Fig. [3J we show 
that the width of the cluster of data points from Protocol 
2 narrows with increasing system size, but the approxi- 
mate linear relationship between the average zj and </>j is 
maintained. Images of five representative packings from 
the scatter plot in Fig. [2] are displayed in Fig. [3] 

Compositional Order We now describe measurements 
of the compositional and positional order for static pack- 
ings. For the compositional order, we quantify the 
fraction of overlapping pairs (ry < cry) that involve 
two small f ss or large fu particles. A scatter plot 
of fu and f ss versus </>,/ for static packings generated 
from both protocols is shown in Fig. 0] The packings 
from Protocol 1 show no signs of phase separation with 
fss + fu ~ fsi ~ 0.5 for all packings. In contrast, Pro- 
tocol 2 generates static packings with a range of compo- 
sitional order as shown in Fig. [3] (c)-(e). For example, 
at the largest <fij, the system displays macrophase sep- 
aration with f ss + fu ~ 1 and f s i ~ 0. We find simi- 
lar results when we define contacting pairs as those with 
Tij < r m i n o~ij, where r m ; n is set by the first minimum in 
9(r)- 

Bond Orientational Order To quantify positional or- 
der, we calculate the bond orientational order parame- 
ter tpQ, which measures the hexagonal registry of near- 
est neighbors [30|. ipe can be calculated 'locally', which 
does not consider phase information, or 'globally', which 
allows phase cancellations. A polycrystal will yield a rel- 
atively large value for the local bond orientational order 



0.8 



0.6 



_^>0.4 



0.2 



(a) 




o 


' 




o 


■ 




°o 
o 






O o ° 

n o 

o ° 8 

o o ° <P 


' 




° ° o ° ° 


• 


<^&- 


° o o %o„_ °„ <P ° 






1 e 8 * ° o 


o 9m 




%o.o 





0.83 0.84 0.85 0.86 0.87 




0.83 



0.84 



0.85 



0.86 



0.87 



f 



FIG. 5: Scatter plot of the (a) global and (b) local bond 
orientational order parameters, ipg and ipg, versus packing 
fraction for static packings from protocol 1 (squares) and 2 
(circles). 



parameter ipg, even though the global order parameter 
ipg ~ 1/y/Nd, where Nd is the number of polycrystalline 
domains. Eqs. ((4|) (global) and (0 (local) provide ex- 
pressions for the bond orientational order parameters in 
2D. 
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where 6ij is the angle between a central particle i and 
neighbors j and n^ denotes the number of nearest neigh- 
bors of i. Two particles are deemed nearest neighbors if 
their center-to-center separation nj < r m [ n aij. 

The results for the global and local bond orientational 
parameters ipg and ip\ are shown in Fig. [5j The static 
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FIG. 8: Eigenvectors corresponding to the modes with fre- 
quencies near the (a) first and (b) second peaks in the den- 
sity of states D(uj) for monodisperse packings with zj — 6 
and 4>j ~ <^ x tai for N = 256. The size of the eigenvector com- 
ponent for each particle is proportional to the length of the 
vector associated with each particle. 



FIG. 6: Density D(uS) of normal mode frequencies uj for 
N = 1024 bidisperse frictionless disk packings obtained us- 
ing Protocols 1 and 2 as a function of the contact number 
at jamming onset for zj ~ 4.0 (black), 4.0 < zj < 4.1 
(red), 4.1 < zj < 4.2 (green), 4.3 < zj < 4.4 (blue), and 
4.5 < zj < 4.6 (violet). The inset shows the same data ex- 
cept that it focuses on low frequencies uj < 1 and includes 
power-law fits to D(uj) ~ uj a as dashed lines. 




FIG. 7: Density D(uj) of normal mode frequencies uj for N = 
1024 monodisperse frictionless disk packings obtained using 
Protocol 1 as a function of the contact number at jamming 
onset for 4.1 < zj < 4.2 (green), 4.5 < Zj < 4.6 (violet), 
4.9 < zj < 5.0 (cyan), 5.4 < zj < 5.5 (magenta), and zj ~ 
6.0 (orange). The inset shows the same data except that it 
focuses on low frequencies uj < 1 and includes power-law fits 
to D(uj) ~ ui a as dashed lines. 



packings obtained from Protocol 1 possess only local 
bond orientational order with ip l 6 ss 0.55 as found in dense 
liquids [3fJ, and ipQ ~ 1/y/N. Further, there is no cor- 
relation between the packing fraction </>j and global or 
local bond orientational order. In contrast, for the phase- 
separated and partially crystalline packings from Proto- 
col 2, we find that there is a strong positive correlation 
between ip l e and cj)j and a somewhat weaker correlation 
between ipg and <j>j. 

The static packings from Protocols 1 and 2 have differ- 
ent structural properties. Those from 1 are amorphous 
and possess similar structural properties even though 
they exist over a range of packing fraction. In con- 
trast, there is a positive correlation between composi- 
tional and positional order and packing fraction for the 
phase-separated and partially crystalline packings from 
Protocol 2. We will now describe the mechanical prop- 
erties of the static packings including the spectrum of 
normal modes and static shear modulus as a function of 
contact number and order. 

Spectrum of Normal Modes The spectrum of nor- 
mal modes provides significant insight into the structural 
and mechanical properties of mechanically stable pack- 
ings [13| • For example, there is evidence that the low- 
frequency region of the spectrum controls the static shear 
response of jammed packings [31|. To calculate the spec- 
trum, we diagonalize the dynamical matrix of all possible 
second derivatives with respect to particle positions eval- 
uated at positions of the static packing — assuming that 
no existing contacts break and no new contacts form [32| ■ 
This yields 2N' — 2 nontrivial eigenvalues e, after ac- 
counting for translational invariance. We consider here 
only mechanically stable packings, and thus all 2N' — 2 
of the eigenvalues are nonzero [33[ . 

The density D(lj) of normal mode frequencies Ui = 
yJei/N, or density of states (DOS), is given by D(uj) — 
(JV(w + Su) — N(id))/6u), where N(u>) is the number of 
modes with frequency less than or equal to uj. The den- 
sity of states D(oj) for packings of bidisperse frictionless 




FIG. 9: Power-law exponent a for the scaling of the density 
of states with frequency in the limit uj — > (D{ui) ~ ui a ) as a 
function of contact number at jamming onset zj for bidisperse 
(circles) and monodisperse (squares) packings. (The error 
bars indicate the error in a from least-squares analysis.) The 
dashed line is a fit to Eq.[7J(with a — 0.17), which interpolates 
the data between the limiting values a = at zj = Zi so = 4 
and q = 1 (Debye behavior) at zj — z xta i = 6. The solid line 
is Eq. [7] with o = 0. 
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FIG. 10: Static shear modulus G versus the deviation in 
packing fraction from the jamming onset A0 = <f> — <j>.j for 
static packings at (zj) = 4.0 (circles), 4.15 (diamonds), 4.35 
(left triangles), and 4.55 (right triangles). The long dashed 
(dot-dashed) line has slope 0.5 (0.4). The inset shows the 
power-law scaling exponent j3 for the static shear modulus 
(G ~ (A(f))P) versus the contact number zj at jamming. 



disks is shown in Fig. [6] as a function of the contact num- 
ber at jamming onset zj. As in previous studies [17] . we 
find that for isostatic systems with zj ~ 4, D(uj) pos- 
sesses a nearly constant regime at low frequencies, which 
signals an abundance of low-frequency modes compared 
to ideal Debye behavior (where D(tu) ~ u as ui — > 0) 
for ideal 2D harmonic solids. For the micro- and macro- 
phase separated bidisperse packings generated using Pro- 
tocol 2 with zj > 4.1, the density of states develops 
two other interesting features. First, D(u>) develops two 
strong peaks near ui ~ 1.0 and 1.6 instead of a single 
broad peak centered near w w 1.4 for isostatic amor- 
phous systems. (We will see below that these peaks are 
associated with crystallization.) Second, we observe that 
as zj increases and the packings become hyperstatic, the 
weight in D(ui) at low frequency (w < 0.3) decreases. As 
shown in the inset to Fig. [51 the density of states scales 
as a power-law 



£>(«) 



(6) 



in the limit uj — > with a scaling exponent a that 
varies continuously with contact number zj as shown 
in Fig. [5] (See Appendix [B] for a discussion of the 
system-size dependence of the exponent a.) Note, how- 
ever, that the plateau in the density of states remains 
largely unchanged in the intermediate frequency regime 
0.3 < uj < 1 over a wide range of zj, which implies that 
some of the remarkable features of jamming in isostatic 
systems also hold for hyperstatic systems. 



To test the generality of the results for the density of 
states, we also calculated D(ui) for monodisperse friction- 
less disk packings generated using Protocol 1 as shown 
in Fig. [71 The density of states for monodisperse systems 
displays similar features to that for bidisperse systems. 1. 
A plateau in D(uj) exists at low to intermediate frequen- 
cies for nearly isostatic systems. 2. Strong distinct peaks 
are located near u — 1.4 and 2.25 for hyperstatic pack- 
ings. Eigenvectors that correspond to the two peak fre- 
quencies are visualized in Fig. [HI 3. A power-law regime 
D(ui) ~ u) a develops in the uj — > limit for hyperstatic 
packings. The exponent a varies continuously with zj 
with a similar functional dependence to that for bidis- 
perse systems as shown in Fig. [9l A notable difference 
between bidisperse and monodisperse systems is that a 
continuous power-law regime in D(uS) persists to higher 
frequencies (oj ~ 1) for monodisperse compared to bidis- 
perse systems. 

The dependence of the scaling exponent a on zj is 
displayed for all bidisperse and monodisperse packings 
(binned by zj) in Fig. [9] We find that a increases mono- 
tonically with Zj and use the suggestive empirical form 



(d-1) 



zj 



Zxtal 



+ a(zj -Ziso)(zj- Zxtal), (7) 



where a is a fitting parameter, to describe the data be- 
tween the limiting values a = at z j = z iso and a — d — 1 
(Debye behavior) at zj = z x tai- The continuous increase 
in a from to 1 as the contact number increases sug- 
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FIG. 11: The contact number z.j as a function of a, where 
the condition rij < (1 + a)erij determines whether particles 
i and j are in contact. The packings shown are TV = 1024, 
4>j = 0.837 (circles); N = 1014, (j>j = 0.892 (squares); and 
N = 2390, 4>j = 0.897 (diamonds). 
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FIG. 12: Contact number zj versus packing fraction <f>j for 
the same data in Fig. [2] and an additional set of packings 
obtained from thermalizing the configurations in Fig. [2] with 
4>.j > 0.86 and then identifying the nearest packing. The 
variation in zj increases with d>j. 



gests a different scenario for the behavior of the jamming 
transition as a function of zj and positional order com- 
pared to the first-order-like transition found as the sys- 
tem compacts above random clo se p acking in simulations 
of frictional granular materials [25| . 

Static Shear Modulus To measure the static linear 
shear modulus G, we slightly deform the system by ap- 
plying an infinitesimal simple shear strain 7 (along the 
x-direction with gradient in the y-direction) , allowing the 
system to relax via energy minimization at fixed strain, 
and then measuring the resulting shear stress response, 
G = dExy/d/y. In Fig. [TQl we show the shear modulus 
versus the amount of compression A0 = (f>— <pj for bidis- 
perse packings obtained from Protocols 1 and 2 at several 
values of zj. We find generally that in the limit A(f> —> 
the static shear modulus scales as a power-law with A</>: 



G = Go(A^)' 3 , 



(8) 



where the scaling exponent j3 (and prefactor Go) depend 
on zj. As shown in Fig.QIJl j3 decreases steadily from 0.5 
to 0.4 as the contact number zj at jamming increases. 
Note that j3 = 0.5 for zj = zi S o was obtained in previous 
work on isostatic packings [17j. The results in Fig. [TU1 
suggest that the critical behavior (e.g. power-law scaling 
of the shear modulus) found in jammed isostatic systems 
persists when the jamming onset is hyperstatic. Further 
studies are required to determine whether the scaling ex- 
ponent for the static shear modulus can be varied over 
the full range from 0.5 to 0. 



IV. CONCLUSIONS 

Using computer simulations, we generated a large li- 
brary of mechanically stable packings of bidisperse, fric- 



tionless disks that span a wide range of contact number 
from zj — Ziso — 4 to z x tai = 6 and packing fraction at 
jamming from <pj ~ 0.84 to near </>xtai- We find that there 
is an amorphous, isostatic branch of packings that spans 
a finite range in packing fraction in the large-system 
limit. Over this range of packing fraction, these packings 
are amorphous with no correlation between bond orien- 
tational order or compositional order and (f>j. We also 
find a branch of phase-separated and partially crystalline 
packings for which the compositional and positional or- 
der increase with <pj. In addition, we characterize the 
mechanical properties of the static packings by measur- 
ing the spectrum of normal modes and the static shear 
modulus. We find that the mechanical properties of the 
packings vary continuously as the contact number and 
structural and compositional order at jamming onset in- 
crease from their isostatic values. In particular, we find 
that the static shear modulus scales as a power-law in the 
amount of compression, G ~ (A0)°, and that the low- 
frequency density of states scales as a power-law in fre- 
quency, D(oj) ~ uj a , and both a and /3 vary continuously 
with contact number at jamming onset. These findings 
emphasize that jamming behavior in systems with purely 
repulsive contact potentials occurs over a range of con- 
tact numbers, not just near zj = Z[ so [34H36l | . In future 
studies, we will investigate the relationship between the 
scaling exponents a and /3, which is likely an important 
feature of jamming in hyperstatic systems. 
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Appendix A: Error analysis of contact number 

In this appendix, we study how sensitive the contact 
number z j is to the definition of whether two particles are 
in contact. In Fig. II H we show zj versus log 10 a where 
two disks i and j are considered in contact (or overlap- 
ping) if Tij < (1 + a)<Jij for three representative config- 
urations: N = 1024, (J) j = 0.837 (circles); N = 1014, 
cf)j = 0.892 (squares); and N = 2390, <pj = 0.897 (dia- 
monds). We see that the contact number is well-defined 
for amorphous configurations at low packing fractions, 
i.e. the contact number is constant over a wide range of 
a that determines whether two particles are in contact. In 
contrast, for packings with large </>j and significant order 
as shown in Fig.[3](e), the contact number varies contin- 
uously with a down to the numerical precision of the par- 
ticle positions in the simulations (a m - m ~ 10~ 8 ). Thus, 
at the current numerical precision of the simulations it 
is difficult to determine zj accurately for the partially 
ordered and ordered configurations. To test the robust- 
ness of the contact numbers, we also added weak thermal 
fluctuations to the packings with 4>j > 0.855 in Fig.[5]for 
times significantly shorter than the structural relaxation 
time, and then found the nearest static packing. This 
data, shown by the small filled symbols in Fig. [T^l pos- 
sess surprisingly small contact numbers and begin to fill 
in the region at large <pj and small zj. As a result, we 
only include configurations in Fig.[5]that possess plateaus 
in zj versus a over a range a m i n < a < a max of at least 
two orders of magnitude. 



Appendix B: Robustness of the Density of States 

In this appendix, we test the robustness of our mea- 
surements of the the density of states D(w) by (1) study- 
ing the system-size dependence of the accumulated fre- 
quency distribution N(uS) and (2) comparing D(uj) for 
hyperstatic packings at jamming onset with contact num- 
ber zj to that for overcompressed packings at the same 
contact number z = Zj. 

To eliminate noise from numerical differentiation, 
we calculate the accumulated distribution N(uj) = 



J Q D(ui')duj' (number of modes with frequency less than 
or equal to u>). For reference, we first show N(u) for 
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FIG. 13: Number N(uj) of normal modes of the dynamical ma- 
trix with frequency less than or equal to u) for monodisperse 
packings at jamming onset with Zj ~ 6 and (/>j — 0xtai and 
N = 16 (circles), 64 (squares), 256 (diamonds), 1024 (upward 
triangles), 2304 (leftward triangles), and 6400 (downward tri- 
angles). The solid line has slope 2. 



monodisperse packings at jamming onset with zj ~ 6 
and 4>.J — </>xtai as a function of system size for N — 16 to 
6400. The crystalline systems show robust Debye power- 



law scaling N(u>) 



at low frequency for all system 



sizes. N(u>) for bidisperse packings at jamming onset is 
shown in Fig. [TJ] for 4.4 < Zj < 4.5 as a function of 
system size. iV(w) displays a power-law scaling with an 
exponent that approaches 1 + a = 1.16 > 1 in the large- 
system limit. Similar robust scaling exponents are found 
for all z j. 

Distinctive features of the density of states D(ui) for 
hyperstatic bidisperse packings at jamming onset are the 
power-law scaling of D(cj) ~ uj° at the lowest frequencies, 
where a varies continuously with zj, and the persistence 
of the plateau in D(lo) at intermediate frequencies over a 
range of zj. Do highly compressed packings display these 
same features? In Fig. 1151 we compare D(u>) for hyper- 
static packings at jamming onset with 4.4 < zj < 4.5 
and overcompressed packings in the same range of con- 
tact number z ~ zj. For the overcompressed packings, 
we find that D(lo) ~ u a , with a — 1, while a « 0.16 
at the lowest frequencies with a crossover to a plateau 
at intermediate frequencies for the hyperstatic packings 
at jamming onset. Thus, hyperstatic packings at jam- 
ming onset possess significantly more low-frequency nor- 
mal modes than overcompressed systems at the same con- 
tact number as shown in the inset to Fig. 1151 
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FIG. 14: Number N(ui) of normal modes of the dynamical 
matrix with frequency less than or equal to u) for bidisperse 
packings at jamming onset generated using Protocol 2 with 
4.4 < zj < 4.5 and N = 512 (circles), 1024 (squares), 2048 
(diamonds), and 4096 (triangles). The solid (dashed) line has 
slope 1.16 (1). 
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FIG. 15: The density of normal modes D(cj) with frequency 
u) for bidisperse packings at jamming onset generated using 
Protocol 2 with 4.4 < zj < 4.5 (blue line) and overcompressed 
packings with contact number z in the same range (red line) . 
The dashed lines in the inset have slope 0.16 and 1. 
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